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Abstract.  Statistical  Model  Checking  (SMC)  is  a  technique,  based  on 
Monte-Carlo  simulations,  for  computing  the  bounded  probability  that  a 
specific  event  occurs  during  a  stochastic  system’s  execution.  Estimating 
the  probability  of  a  “rare”  event  accurately  with  SMC  requires  many 
simulations.  To  this  end,  Importance  Sampling  (IS)  is  used  to  reduce  the 
simulation  effort.  Commonly,  IS  involves  “tilting”  the  parameters  of  the 
original  input  distribution,  which  is  ineffective  if  the  set  of  inputs  causing 
the  event  (i.e.,  input-event  region)  is  disjoint.  In  this  paper,  we  propose  a 
technique  called  Semantic  Importance  Sampling  (SIS)  to  addresses  this 
challenge.  Using  an  SMT  solver,  SIS  recursively  constructs  an  abstract 
indicator  function  that  over-approximates  the  input-event  region,  and 
then  uses  this  abstract  indicator  function  to  perform  SMC  with  IS.  By 
using  abstraction  and  SMT  solving,  SIS  thus  exposes  a  new  connection 
between  the  verification  of  non-deterministic  and  stochastic  systems.  We 
also  propose  two  optimizations  that  reduce  the  SMT  solving  cost  of  SIS 
significantly.  Finally,  we  implement  SIS  and  validate  it  on  several  prob¬ 
lems.  Our  results  indicate  that  SIS  reduces  simulation  effort  by  multiple 
orders  of  maganitude  even  in  systems  with  disjoint  input-event  regions. 


1  Introduction 

Many  systems  deployed  in  the  real-world  are  stochastic,  i.e.,  their  behavior  de¬ 
pends  on  random  inputs  (e.g.,  sensor  readings,  task  execution  times,  etc.)  As 
these  systems  become  more  complex,  there  is  a  growing  demand  for  efficient  and 
precise  techniques  to  verify  correctness  of  their  behavior.  In  this  paper,  we  target 
a  common  verification  problem  -  estimating  the  probability  of  an  event  r  (e.g., 
some  sort  of  failure)  during  the  execution  of  a  stochastic  system  Af.  Analytic 
solutions  to  this  problem  (e.g.,  probabilistic  model  checking,  see  Section  2)  do 
not  scale  to  many  real-world  systems  due  to  complexity.  We  focus  on  an  al¬ 
ternate  approach  called  Statistical  Model  Checking  (SMC)  [11],  which  relies  on 
Monte-Carlo-based  simulations  to  solve  this  verification  task  more  scalably. 

*  This  material  is  based  upon  work  funded  and  supported  by  the  Department  of  De¬ 
fense  under  Contract  No.  FA8721-05-C-0003  with  Carnegie  Mellon  University  for  the 
operation  of  the  Software  Engineering  Institute,  a  federally  funded  research  and  de¬ 
velopment  center.  This  material  has  been  approved  for  public  release  and  unlimited 
distribution.  DM-0001772 
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(a)  (i)  (c) 


Fig.  1.  Example  of  SIS;  f(x)  =  original  input  distribution;  g(x)  =  tilted  distribution; 
g*(x )  =  distribution  produced  by  SIS. 


SMC  produces  two  results  -  the  estimate  p  of  the  probability  p  of  r  and  a 
measure  e  of  precision  of  p.  The  key  challenge  in  SMC  is  “simulation  explosion” 

-  the  number  of  simulations  required  to  achieve  a  high  e  becomes  prohibitively 
large  if  p  is  small  (i.e.,  r  is  rare).  Importance  Sampling  [9]  (IS)  has  been  shown 
to  address  this  challenge.  Suppose  the  random  input  x  to  M  has  distribution 
f[x).  In  IS,  we  first  perform  SMC  under  a  different  input  distribution  g{x)  that 
makes  r  more  likely  (i.e.,  increases  p),  and  then  adjust  the  result  back  to  f(x). 

Traditionally,  importance  Sampling  is  implemented  by  “tiliting”  the  param¬ 
eters  of  the  input  distributions  to  increase  the  likelihood  of  r.  However,  tilting  is 
less  effective  if  the  set  of  inputs  that  cause  t,  i.e.,  the  input-event  region  denoted 
xT ,  is  disjoint.  For  instance,  this  happens  when  analysing  a  program  where  r 
only  occurs  if  the  execution  follows  one  of  several  control-flow  paths,  each  trig¬ 
gered  by  a  distinct  input  range.  Figure  1(a)  shows  such  a  case.  The  actual  input 
distrbution  f(x)  is  uniform  in  the  range  [0, 10],  and  xT  =  [2.99,  3.01] U  [6.99,  7.01]. 
Figure  1(b)  shows  a  titled  distribution  g(x)  uniform  in  the  range  [2,10].  While 
g{x)  make  r  more  likely  than  f(x),  it  still  assigns  positive  weight  to  large  parts 

-  e.g.,  (3.01,  6.99)  -  of  the  input  space  that  do  not  belong  to  xT. 

In  this  paper,  we  address  this  challenge,  and  make  three  specific  contribu¬ 
tions.  First,  we  develop  a  new  technique  to  construct  more  precise  input  distr- 
butions  for  IS  -  such  as  g*(x)  shown  in  Figure  1(c)  -  even  when  the  input-event 
region  is  disjoint.  This  technique,  which  we  call  Semantic  Importance  Sampling 
(SIS),  takes  as  input  a  description  of  M.  and  f(x),  and  recursively  computes  a 
precise  “over-approximation”  of  xT  in  the  form  of  an  abstract  indicator  func¬ 
tion  (AIF).  In  each  step  of  the  recursion,  SIS  constructs  a  verification  condition 
using  A4  and  f(x)  and  checks  its  satisfiability  with  a  SMT  solver  to  eliminate 
parts  of  the  input  space  that  are  not  in  xT.  The  algorithm  outputs  an  AIF  rep¬ 
resented  by  a  set  of  “input  cubes”,  i.e.,  a  disjunction  of  intervals  [?]  over  the 
input  variables  of  M.  Subsequently,  SIS  uses  the  AIF  to  construct  a  precise 
input  distribution,  and  perform  SMC  with  IS.  By  using  the  semantics  of  A4, 
SIS  succesfully  applies  concepts  and  techniques  used  widely  in  the  verification 
of  non-deterministic  systems  (such  as  abstraction,  SMT  solving,  and  verification 
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conditions)  to  the  analysis  of  stochastic  systems.  In  this  way,  SIS  builds  new 
bridges  between  these  two  disciplines. 

The  most  expensive  component  of  SIS  are  the  calls  to  the  SMT  solver.  Our 
second  contribution  is  two  optimizations  to  SIS  that  reduce  the  number  of  SMT 
calls  while  maintaining  correctness.  Finally,  we  implement  SIS  in  a  tool  called 
OSMOSIS  and  use  it  to  verify  a  number  of  stochastic  systems  with  rare  events. 
Our  results  indicate  that  SIS  reduces  the  number  of  simulations  significantly, 
in  some  cases  by  a  factor  of  over  600,  and  verification  time  by  over  two  orders 
of  magnitude.  Furthermore,  our  optimizations  reduce  both  the  number  of  SMT 
calls  and  overall  SMT  solving  time,  typically  by  a  factor  of  2.  All  our  tools  and 
examples  are  available  at  andrew.cmu.edu/~schaki/misc/osmosis.zip. 

The  rest  of  the  paper  is  organized  as  follows.  Section  2  surveys  related  work. 
Section  3  presents  background  definitions  and  concepts.  Section  4  presents  SIS, 
and  Section  5  presents  our  tool  OSMOSIS.  In  Section  6,  we  present  our  experiments 
and  results,  and  in  Section  7,  we  conclude. 


2  Related  Work 

SC:  This  needs  a  lot  of  work.  Like  statistical  model  checking,  probabilistic 
model  checking  [10]  is  also  an  automated,  algorithmic  approach  for  computing 
numerical  properties  of  stochastic  systems.  However,  in  this  approach,  the  sys¬ 
tem  is  modeled  as  a  finite  state  probabilistic  automaton,  e.g.,  a  discrete  time 
Markov  chain  (DTMC),  a  continuous  time  Markov  chain  (CTMC),  or  a  Markov 
decision  process  (MDP)  which  is  exhaustively  explored  in  the  analysis.  The  prop¬ 
erty  is  expressed  as  formula  in  a  temporal  logic,  e.g.,  probabilistic  Computation 
Tree  Logic  (PCTL)  [5].  For  example,  the  system  could  be  a  DTMC  modeling 
the  repeated  throwing  of  a  biased  coin  that  comes  up  “heads”  with  probability 
0.55  and  “tails”  with  probability  0.45.  The  property  could  be  the  “probability 
of  seeing  3  heads  followed  by  4  tails  followed  by  5  heads”.  Probabilistic  model 
checking  then  involves  an  exhaustive  exploration  of  the  DTMC’s  statespace  to 
compute  the  property’s  value.  Typically,  this  involves  the  following  steps:  (i) 
compute  the  set  of  reachable  discrete  states  S  of  an  automaton  constructed  by 
composing  the  system  with  the  property;  (ii)  construct  a  set  of  equations  Q 
whose  solution  corresponds  to  the  steady-state  probabilities  of  S;  and  (iii)  solve 
Q  numerically  and  extract  the  value  of  the  property  from  the  solution.  Proba¬ 
bilistic  model  checking  is  an  active  area  of  research,  involving  both  theoretical 
advancements  [3]  and  practical  tool  development  [7].  It  has  been  used  to  verify 
systems  ranging  from  pacemakers  [1],  root  contention  protocols  [8]  and  biologi¬ 
cal  pathways  [6] .  For  our  approach,  we  used  statistical  model  checking,  since  the 
systems  we  want  to  verify  are  too  complex  to  be  modeled  as  Markov  chains. 


3  Background 

Consider  a  system  A4  with  finite  vector  of  random  inputs  x.  Assume  that  A4  is 
deterministic,  i.e.,  its  behavior  is  fixed  for  a  fixed  value  of  x.  The  SMC  problem 
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P  =  J  lM\=*(x)f(x)dx  P  ^T,^=lIM\=AXi) 


RE(p)  = 


1  -P 
pN 


N  = 


RE(p)  = 
1  ~P 


\J  Var(p) 

m 


P  =  J  /.vi:  .„(x)  ^  ^g(T)d., 


pRE2(p) 

f  IM^#(x)W(x)g(x)dx 


Fig.  2.  Background  definitions.  Var  =  variance;  E  =  expected  value;  W(x)  =  ^2. 


is  to  estimate  the  probability  that  M  satisfies  a  property  <P,  denoted  M  \=  <P, 
given  a  joint  probability  distribution  f{x)  on  x,  i.e.,  to  estimate  p  =  Pr[  M.  |=  <P\. 
We  assume  that  whether  M  \=  <P  under  input  x  can  be  determined  by  simulating 
Af  for  finite  time.  Specifically,  we  assume  that  Af  is  a  program  that  terminates 
under  all  inputs,  and  Af  |=  under  input  x  iff  the  execution  of  Af  under  input 
x  violates  an  assertion  (representing  a  desired  safety  property)  in  A4. 

Let  us  write  x~  f(x)  to  mean  x  is  distributed  by  /( x).  SMC  involves  a  series 
of  Bernoulli  trials,  modeling  each  trial  as  a  Bernoulli  random  variable  having 
value  1  with  probability  p ,  and  0  with  probability  1  —p.  For  each  trial,  a  random 
vector  x~  f(x)  is  generated.  The  system  Af  is  simulated  with  input  x  to  generate 
a  trace  a.  The  trial’s  outcome  is  0  if  $  holds  on  a,  and  0  otherwise. 

Define  an  indicator  function  Im\=&(x)  that  returns  1  if  M.  |=  <P  under  input 
x,  and  0  otherwise.  Then  the  probability  p  can  be  calculated  as  in  Figure  2.  It 
can  be  estimated  as  p  shown  in  Figure  2,  where  N  is  the  number  of  trials,  and 
Xj~/(x).  The  precision  of  p  is  quantified  by  its  “relative  error”  RE{p),  defined 
in  Figure  2.  The  second  row  in  Figure  2  shows  a  known  relationship  [2]  beween 
RE(p),  p  and  N,  and  its  rearrangement  with  N  on  the  left  of  the  equality. 

Importance  Sampling.  Figure  2  shows  that  the  number  of  simulations 
needed  to  achieve  a  fixed  precision  with  SMC  increases  rapidly  as  the  target 
event  becomes  rarer.  Importance  Sampling  [9]  (IS)  has  been  applied  [2]  to  address 
this  challenge  effectively  by  reducing  Var(p).  The  key  idea  behind  IS  is  to  first 
simulate  A4  under  a  different  input  distribution  g(x)  with  lower  variance,  and 
then  mathematically  adjust  the  result  back  to  the  original  distribution  /( x). 
The  third  row  of  Figure  2  shows  the  definition  of  p  using  g{x)  and  the  “weight 
function”  W{x)  =  The  estimator  for  this  form  is: 


N 

P  =  '52IM\=<p(Xi)W{xi)  (1) 

i-1 

where  the  Xi~  g( x).  The  biggest  challenge  in  applying  IS  effectively  is  choosing  a 
“good”  g(x)  that  will  reduce  Var(p).  Typically  this  is  done  by  “tilting”  f(x)  by 
changing  its  distribution  parameters  (mean,  variance  etc.)  However,  as  discussed, 
tilting  is  not  effective  if  is  disjoint  in  the  input  space.  Our  main  contribution, 
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SIS,  is  a  technique  for  constructing  a  good  g(x)  even  in  such  cases.  We  describe 
this  technique  in  detail  in  the  next  section. 


4  Semantic  Importance  Sampling 


To  explain  SIS,  we  begin  with  a  known  result  [2]  that  there  always  exists  an 
optimal  IS  distribution: 


V 


(2) 


for  which  Var(p)  =  0,  i.e. ,  if  IS  is  done  with  g(x)  =  g<>{x),  then  a  single  sample  is 
sufficient  to  compute  p.  However,  there  are  two  challenges  to  using  g°(x)  for  IS: 
(i)  50  (x)  depends  on  p,  the  answer  we  are  actually  looking  for;  and  (ii)  g<>{x)  also 
depends  on  the  indicator  function  Im\=& (x)>  but  since  this  function  represents 
A4  itself,  it  may  be  too  complex  to  represent  analytically. 

The  key  insight  behind  SIS  is  to  construct  an  abstract  indicator  function 
(AIF)  such  that:  (i)  IM\=$(X)  =  1  =>  =  b  and  (U) 

is  simple  enough  to  represent  analytically.  Note  that  the  AIF  represents 

an  over-approximation  of  the  set  of  inputs  under  which  A4  satisfies  <P.  This  AIF 
induces  the  following  IS  distribution  and  weight  function: 


TM)=*(x)f(x) 

p* 

=  fix)  =  f(x)p*  =  p* 
9*{x)  IM^{x)f{x) 


(3) 

(4) 


where  p*  =  E[I^^(x)  =  1]  is  the  probability  that  for  an  input  x  ~  /( x), 
Im\=$(.x)  =  4-  Note  that  as  I*M ,  ^(a;)  approaches  Im\=<p(x)i  9*(x )  also  aP" 
proaches  g<>{x).  In  the  limit,  =  Im\=<p{x)  implies  g*(x)  =  g<>{x). 

Probability  Estimation  and  Relative  Error  in  SIS.  Substituting  W*(x) 
from  (4)  into  (1),  we  get  the  SIS  estimator: 


1  .  .  1  ,  p* 

P=  „'52IM\=*(xi)W*{xi)  =  - At 

JV  i= l  JV  i= i 

where  Xi~g*{x).  Note  that,  from  (3),  if  Xi~g*(x),  then 
p*  is  a  constant.  Hence,  (5)  simplifies  to: 

*  N 

p  =  —  ^  lM\=&{xi )  =  P*  x  praw,  where 


i—1 


Praw  — 


1 

N 


N 

'y  '  lM\=®{.xi) 


i=  1 


(5) 

=  1.  Also, 

(6) 
(7) 
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is  the  estimated  probability  that  for  an  input  Xi~g*{x),  lM\=&(xi)  =  1)  he.,  x%  is 
an  input  where  A4  \=  <P  holds.  Since  praw  is  an  unweighted  average  of  Bernoulli 
random  variables,  its  relative  error  can  be  estimated  [2]  as: 


RE(piaw) 


1 

VPi'&wN 


(8) 


Furthermore,  from  (6),  since  p  =  p*  x  praw,  and  p*  is  a  constant,  the  relative 
error  for  p  is  the  same  as  the  relative  error  of  praw,  he-,  RE{p)  =  i?if(praw)- 


4.1  The  SIS  Algorithm 

The  SIS  algorithm  involves  the  following  steps: 

1.  Use  iterative  abstraction-refinement  to  construct  the  AIF 

2.  Calcuate  p*  and  g*(x). 

3.  Use  SMC  to  estimate  praw  with  desired  RE(p)  =  RE(pi-aw),  using  g*(x)  to 
draw  random  inputs.  Output  p  =  p*  x  praw- 

The  core  of  SIS  is  Step  1,  the  generation  of  the  AIF.  We  describe  this  in 
the  following  sections  by  first  discussing  our  representation  of  the  AIF,  then 
describing  the  abstraction-refinement  algorithm. 

AIF  as  a  Cube  Set.  We  assume  that  the  input  x  to  A4  is  a  vector  of  M 
independent1,  but  not  necessarily  identically  distributed  random  variables.  For 
each  dimension  Xi  in  x,  Ei(x)  is  the  Cumulative  Distribution  Function  (CDF), 
F“1(u)  is  the  inverse  CDF  (or  quantile  function),  and  Ui  =  F~1(xi)  is  the 
quantile  domain  variable.  Let  £  be  a  M-dimensional  axis-aligned  input  domain 
hypercube  defining  an  interval  [li,  hi]  on  each  input  variable  Xi  for  1  <  *  <  M .  We 
also  define  the  quantile  domain  hypercube  c  defined  by  the  ranges  [F)(Ij),  F)(/ij)] 
for  each  dimension.  We  use  the  notation  c  =  F(^)  and  £  =  F-1(c)  to  transform 
cubes  between  input  and  quantile  space.  We  will  use  the  terms  input  cube  and 
quantile  cube  to  refer  to  cubes  in  the  input  and  quantile  spaces,  respectively. 
When  the  term  cube  is  used  without  qualification  we  will  assume  quantile  cubes. 
We  can  now  represent  the  AIF  in  terms  of  a  qualtile  cube  set  C*  as: 


T *  /  n  _  f  1  if  3c  G  C*  |  F(x)  €  c 

|  0  otherwise 


(9) 


where  Vc  €  C* ,  Im\=®{x)  =  1  4-  i  £  c. 

Cube  Splitting.  Let  £u  be  the  input  cube  defining  the  support  of  the  input 
distribution  function  f{x).  The  corresponding  qunatile  domain  cube  cjj  =  F(l]u) 
will  have  a  range  of  [0, 1]  on  each  dimension.  We  call  this  the  level-0  cube.  We 
write  c/i  to  mean  the  cube  formed  by  spliting  the  interval  on  iq  in  c  in  half,  and 
retaining  only  the  upper  half.  Similarly,  c/i  is  the  result  of  a  similar  operation 
where  the  lower  half  of  the  interval  is  retained.  Note  that  we  can  split  on  the 


1  Non-independent  random  inputs  y  are  replaced  by  a  function  h(x)  of  independent 
random  variables  x,  which  is  folded  into  to  yield  Im\=& (h(x)). 
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(1)  CubeSet  aifGen (SMT  I, Cube  c) 

(2)  { 

(3)  if  (Solved,  ==  UNSAT)  return  0; 

(4)  if  (level  (c)  ==  Lmax)  return  {c}; 

(5)  int  k  =  (level  (c)/G)  "/,  M ; 

(6)  Cube  Co  =  c/k ;  Cube  ci  =  c/k; 

(7)  return  aifGend,  Co)  U  aifGen(/,  Ci)  ; 

(8)  } 

Fig.  3.  Basic  AIF  Generation  Algorithm;  G=variable  grouping  factor,  M=number  of 
input  variables,  I/max=recusion  depth  limit,  Solve  =  satisfiability  check  via  SMT  solver. 


same  variable  multiple  times.  A  level-fc  cube  is  the  result  of  k  splits  on  the  level-0 
cube.  For  example  if  cjj  is  the  level-0  cube,  then  Q7/I/I  is  the  level-2  cube  in 
which  the  interval  for  u\  is  [0.5,0.75].  After  each  split,  the  probability  that  a 
input  drawn  from  f(x)  falls  in  the  result  is  halved.  Thus,  the  probability  of  an 
input  drawn  from  f(x)  falling  in  a  level-A:  cube  is 

Iterative  Abstraction-Refinement.  Generation  of  the  AIF  (x)  is 

performed  recursively  through  the  hierarchical  use  of  an  SMT  solver.  The  basic 
algorithm  aifGen  is  shown  in  Figure  3.  It  takes  as  input  the  SMT  representation 
of  the  indicator  function  /,  and  the  input  cube  c  over  which  to  generate  an 
abstaction.  Constant  Lm ax  is  the  maximum  recursion  depth.  It  returns  the  subset 
of  level- Amax  cubes  in  C*  within  cube  c.  C*  representing  the  AIF  as  defined  in 
(9)  can  then  be  deterimined  by  calling  aifGen,  passing  the  level-0  cube  cjj  as  c. 

The  algorithm  works  as  follows.  At  Line  3,  the  SMT  solver  is  applied  to  the 
model  /  over  the  cube  £  =  i7'_1(c).  The  cube  is  applied  to  the  model  by  modifying 
the  assertions  in  the  model  corresponding  to  the  intervals  on  the  input  variables. 

The  SMT  solver  can  return  SAT,  UNSAT  or  UNKNOWN  (e.g.,  if  it  times 
out).  If  the  result  is  UNSAT,  then  M  \=  <P  does  not  hold  in  the  input  space 
described  by  c,  and  so  it  returns  the  empty  set.  If  the  result  is  SAT  or  UN¬ 
KNOWN,  we  continue  with  the  rest  of  the  algorithm.  While  an  UNKNOWN 
result  will  reduce  the  efficiency  of  the  algorithm,  the  result  will  still  be  sound. 

At  Line  4,  the  level  of  the  current  cube  c  is  checked  against  the  specified 
maximum  recursion  depth  Lmax.  If  we  are  at  that  maximum  recursion  depth,  we 
simply  return  the  set  containting  just  the  cube  c. 

At  Line  5,  we  choose  an  input  variable  index  on  which  to  split  the  cur¬ 
rent  cube.  In  our  current  implementation,  we  simply  cycle  through  the  variables 
round-robin  by  using  the  current  level  modulo  the  total  number  of  input  vari¬ 
ables  M.  Integer  division  by  a  variable  grouping  factor  G  allows  us  to  choose  the 
same  variable  G  levels  in  a  row  before  moving  to  the  next  variable.  It  is  possi¬ 
ble  that  other  methods  of  chosing  the  spliting  order  may  lead  to  more  efficient 
abstractions,  however  we  have  not  yet  explored  this  area. 

At  Lines  6-7,  we  split  the  cube  c  around  the  selected  variable  U}~  forming  the 
cubes  Co,  and  C\  for  the  lower  and  upper  half  of  the  CDF  interval  on  variable  Uk 
in  c.  We  then  recufrsively  call  the  generation  algorithm  on  those  two  sub-cubes 
and  return  the  union  of  the  cube  sets  returned  by  each  call. 
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Calculation  of  p*.  All  cubes  in  the  set  C*  returned  by  aifGen  are  level- 
Lmax  and  non-overlapping.  As  a  result,  p*  is  calculated  from  the  number  of  cubes 
in  the  set,  and  the  recursion  depth  limit  as: 


P 


* 


\C*\ 

2-^max 


(10) 


which  represents  the  probability  that  an  input  from  the  original  distribution 
f(x)  will  fall  in  the  abstract  indicator  function  I*M ^_$(x). 


4.2  Optimized  AIF  Generation 

The  most  expensive  component  of  aifGen  are  the  calls  to  Solve.  We  now  present 
two  optimizations  that  can  reduce  the  number  of  calls. 

Optimization  1:  Skip  on  UNSAT.  Consider  the  algorithm  in  Figure  3. 
Notice  that  at  the  point  where  we  split  the  cube  at  Line  6,  we  already  know 
that  cube  c  is  not  UNSAT.  The  means  that  if  one  of  the  child  cubes  Co  or  ci 
is  UNSAT,  the  other  one  must  be  SAT2.  To  take  advantage  of  this,  we  modify 
the  algorithm  to  take  an  additional  boolean  argument  assumeSAT  indicating  we 
should  skip  the  call  to  Solve  and  assume  it  returns  SAT  when  assumeSAT  is 
true.  Then  we  make  the  first  recursive  call  on  Co  with  assumeSAT  set  to  false.  If 
this  call  returns  the  empty  set,  then  the  result  for  that  half  was  UNSAT,  and 
we  pass  true  for  assumeSAT  when  making  the  recurisve  call  on  ci,  otherwise  we 
make  the  recusrive  call  with  assumeSAT  set  to  false  and  execute  Solve  as  normal. 

Optimization  2:  Counter-Example  Reuse.  A  second  optimization  is 
possible  by  making  use  of  the  counter-example  returned  by  Solve  when  the 
result  is  SAT.  In  this  case,  we  assume  that  Solve  returns,  as  counter-example, 
a  cube  Ui  containing  a  satisfying  solution.  We  convert  to  a  quantile  space 
cube  Cd  =  F(£d)-  If  Cd  is  completely  contained  by  one  of  the  child  cubes  in  the 
recursive  call,  we  can  skip  the  call  to  Solve  for  that  call.  We  require  Cd  to  be 
completely  contained  since  the  counter-example  cube  Ui  returned  by  Solve  is  a 
cube  in  which  there  exists  a  solution  to  the  SMT  formula,  but  not  all  points  in 
the  cube  are  necessarily  a  solution.  In  most  cases  Cd  will  be  contained  by  one 
or  the  other  of  the  child  cubes  in  the  resursive  calls,  but  it  is  possible  that  Cd 
could  fall  on  an  edge  and  thus  not  be  applicable  to  either  recursive  call.  In  this 
case,  it  is  still  possible  that  Optimization  1  can  apply.  We  assume  that  Solve 
will  return  the  empty  cube  0  when  the  result  is  UNKNOWN  which  will  supress 
use  of  this  optimization  for  the  child  invocations. 

Optimized  AIF  Generation  Algorithm.  Figure  4  shows  the  fully  opti¬ 
mized  abstract  indicactor  function  incorporating  both  of  the  optmizations  dis¬ 
cussed  above.  Line  3  tests  for  conditions  that  allow  us  to  skip  the  SMT  check. 
In  the  case  that  we  are  skipping  a  check,  we  can  pass  the  existing  Cd  to  the  child 
recursive  calls  since  it  may  apply  to  one  of  those  calls  as  well.  When  doing  the 
SMT  check  with  Solve  at  Line  4,  we  include  an  additional  return  parameter  Ui 

2  It  could  be  UKNOWN  if  result  from  cube  c  is  UKNOWN,  but  without  loss  of 
soundness  we  treat  an  UNKNOWN  as  SAT  for  the  purpose  of  this  optimization. 
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(1)  CubeSet  aifGen(SMT  /, Cube  c, boolean  assumeSAT,Cube  Cd) 

(2)  { 

(3)  if  (lassumeSAT  &&  Cd  !  =  0  &&  !(cdCc))  { 

(4)  if  (Solved,  F_1(c),  kxid)  ==  UNSAT)  return  0; 

(5)  cd  =  F(£d); 

(6)  } 

(7)  if  (level(c)  ==  Lmax)  return  {c}; 

(8)  int  fc  =  (level  (c)/G)  7  M; 

(9)  Cube  Co  =  c/k;  Cube  ci  =  c/k; 

(10)  CubeSet  So  =  aifGend,  Co ,  false,  Cd)  ; 

(11)  CubeSet  Si  =  0 ; 

(12)  if  (so  ==  0)  Si  =  aifGend,  Ci ,  true,  Cd)  ; 

(13)  else  Si  =  aifGend,  Ci ,  false,  Cd)  ; 

(14)  return  So U Si ; 

(15)  } 

Fig.  4.  Optimized  Abstract  Indicator  Function  (AIF)  Generation  Algorithm; 
G=variable  grouping  factor,  M=number  of  input  variables,  Lmax=recusion  depth  limit. 


in  which  the  counter-example  cube  is  returned.  We  assume  that  the  empty  set  is 
returned  if  the  result  is  not  SAT.  At  Line  5  we  convert  the  input  space  cube  £<* 
to  a  quantile  space  cube  Cd ■  Lines  12  to  13  implement  Optimization  1.  If  so  = 
then  the  result  of  the  test  for  Co  was  UNSAT  and  we  can  assume  that  the  test 
for  Ci  will  be  SAT. 


4.3  Statistical  Model  Checking 

After  generating  the  AIF  /d=<J5(a;),  and  computing  p*  and  g*{x),  the  last  step 
in  SIS  is  the  actual  SMC.  As  peviously  mentioned,  we  draw  samples  from  the 
distribution  g*(x)  as  defined  in  (3),  then  use  (7)  to  estimate  the  raw  probability 
Prow  and  scale  this  by  the  value  p*  calculated  above. 

No.  Of  Samples.  The  number  of  samples  N*  needed  to  estimate  praw  is: 

N*  _  1  ~  Praw  _  1  ~  P/P  /-q-j 

Praw RE^  {j) raw)  Vl P*  RE^  ij^raw') 


From  (6),  we  know  that  RE(p)  =  RE(praw).  Assuming  small  p  and  p*  p,  the 
speedup  due  to  SIS  can  be  estimated  as: 


AT 

A  _  pRE2(p) 

N*  ~  i -p/p* 

p/p*RE2(praw) 


1  —  p  ^  1 
p*  -p  ~  p* 


(12) 


Random  Input  Generation.  To  generate  a  random  input  from  tj"  (x). 
we  recognize  that  this  is  the  equivalent  of  generating  an  input  from  f(x)  and 
accepting  only  those  that  fall  within  We  do  this  by  first  randomly 

selecting  a  cube  c  from  C*  with  uniform  probability  since  each  cube  has  equal 
probability  of  containing  a  sample  drawn  from  f(x).  We  then  choose  a  uniform 
vector  u  €  c  and  use  the  inverse  CDF  to  generate  the  input  vector  as  x  =  F~l  (u). 
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Fig.  5.  Architecture  of  osmosis  Tool. 


5  Osmosis 

We  implemented  SIS  in  a  tool  called  OSMOSIS.  The  input  to  OSMOSIS  is  a  de¬ 
scription  of  M  in  an  anotated  version  of  C,  with  the  target  property  <P  defined  as 
ASSERT  ()  statements.  OSMOSIS  calculates  the  probability  of  an  ASSERT  ()  failure 
via  SIS,  using  dReal[4]  as  the  backend  SMT  solver. 

Osmosis  Architecture.  Figure  5  shows  the  architecture  of  OSMOSIS.  The 
input  model  is  processed  by:  (i)  gcc  to  generate  a  dynamic  executable;  (ii)  a 
syntactic  extractor  which  looks  for  //Odist  declarations  to  determine  the  input 
space  and  distributions;  and  (iii)  a  verification  condition  generator  that  generates 
an  SMT  formula  corresponding  to  the  C  model.  Then  aifGen  (from  Figure  3  or 
Figure  4)  is  used  to  build  the  AIF  |=<j5  (x) .  This  AIF  is  used  to  calculate  p* 
and  g*{x ),  and  in  conjunction  with  the  dynamically  loaded  executable  for  A4  to 
calculate  praw  and  ILE(praw)-  Finally,  p  is  calculated  using  p*  and  praw- 

Osmosis  Input  Format.  Figure  6(a)  shows  an  example  OSMOSIS  input 
model.  The  anotations  at  Lines  5  and  6  indicate  the  input  to  the  model.  Line  5 
defines  a  random  input  named  “a”  with  a  uniform  distribution  between  0  and  5. 
Line  6  defines  a  random  input  named  “b”  with  a  normal  distribution  with  mean 
3,  standard  deviation  1  which  has  been  censored  to  be  between  0  and  5.  Where 
apprortiate,  we  refer  to  the  model  input  collectively  as  he  vector  x. 

There  are  two  special  functions/macros  in  OSMOSIS  models:  (i)  ASSERTO 
defines  a  condition  that  is  expected  to  be  true;  and  (ii)  INPUT_D()  accesses  a 
random  input  declared  in  an  annotation.  The  suffix  _D  on  INPUT_D()  indicates 
the  return  type  of  double.  In  Figure  6(a),  Lines  8  and  9  access  inputs  “a”  and 
“b”  and  place  them  in  C  variables  also  named  “a”  and  “b” .  Some  computations 
are  performed  on  lines  10  and  1,  then  finally  an  assertion  is  made  on  Line  13.  The 
include  on  Line  1,  allows  the  model  include  the  special  functions  to  be  compiled 
by  a  standard  compiler  such  as  gcc  for  use  in  the  SMC  phase. 

SMT  Generation.  In  order  to  implement  Solve,  OSMOSIS  translates  the  C 
model  into  a  verification  condition  represented  as  an  SMT  formula  <p,  which  is  in 
essence,  a  representation  of  the  indicator  function  IM^=<p(x),  i.e.,  any  input  value 
x  satisfies  ip  iff  Im\=<s>{ x)  =  1.  In  constructing  p,  stochastic  inputs  defined  by  the 
@dist  annotations  in  the  C  model  use  the  same  variable  name  as  the  declaration. 
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(1)  #include  "osmosis_model .h" 

(2) 

(3) 

(4)  //@dist  a=uniform(min=0,max=5) 

(5)  //Odist  b=normal(mean=3,std=l , 

min=0,max=5) 

(6)  void  model () 

(7)  { 

(8)  double  a  =  INPUT_D("a") ; 

(9)  double  b  =  INPUT_D("b") ; 

(10)  double  c  =  a  +  b; 

(11)  double  d  =  (a  -  b)/2.0; 

(12) 

(13)  ASSERT(sin(c)*cos(d)  <=0.999); 

(14)  > 


(1)  (set-logic  QF_NRA) 

(2)  (declare-fun  a  ()  Real) 

(3)  (declare-fun  b  ()  Real) 

(4)  (declare-fun  a_l  ()  Real) 

(5)  (declare-fun  b_l  ()  Real) 

(6)  (declare-fun  c_l  ()  Real) 

(7)  (declare-fun  d_l  ()  Real) 

(8)  (assert  (>=  a  0)) 

(9)  (assert  (<=  a  5)) 

(10)  (assert  (>=  b  0)) 

(11)  (assert  (<=  b  5)) 

(12)  (assert  (=  a_l  a)) 

(13)  (assert  (=  b_l  b)) 

(14)  (assert  (=  c_l  (+  a_l  b_l))) 

(15)  (assert  (=  d_l  (/  (-  a_l  b_l)  2))) 

(16)  (assert  (not  (<=  (*  (sin  c_l) 

(cos  d_l) )  0.999))) 

(17)  (check-sat) 

(18)  (exit) 


(a) 


if  (a  >  b) 

a  =  cos(a*b) ; 


(b) 

(assert  (=  _C1  (>  a_l  b_l))) 

(assert  (or  (not  _C1)  (=  a_2  (cos  (*  a_l  b_l))))) 
(assert  (or  _C1  (=  a_2  a_l))) 


(c) 


(d) 


Fig.  6.  (a)  OSMOSIS  Input  Example;  (b)  SMT  for  OSMOSIS  Input  Example;  (c)  a  con¬ 
ditional  statement;  and  (d)  its  translation  to  SMT. 


The  model  is  also  converted  to  single-static-assignment  form  so  that  each  local 
variable  is  assigned  once.  Finally,  ASSERTO  conditions  are  negated  since  we  are 
interested  in  testing  if  there  are  any  inputs  that  can  result  in  an  assertion  failure. 

Assignments  become  equalities.  Conditional  (if)  statements  are  translated 
by  generating  a  variable  for  the  condition,  then  translating  both  branches  as 
consequences  of  implications  of  the  condition,  or  the  compliment  of  the  condi¬ 
tion.  If  there  are  differing  numbers  of  assignments  to  a  variable  in  the  branches, 
then  an  additional  assertion  is  added  to  reconcile  the  generation  numbers  of 
the  variables.  For  example,  the  C  statement  in  Figure  6(c)  generates  the  SMT 
assertions  in  Figure  6(d).  Loop  (while  and  for)  statements  are  unrolled  and 
must  include  an  annotation  to  indicate  the  maximim  loop  count.  Note  that  the 
construction  of  p  is  effective  and  linear  in  the  size  of  the  model. 

Figure  6(b)  shows  the  p  generated  from  the  M  given  in  Figure  6(a).  Line 
8  through  11  define  the  intervals  in  the  stochastic  inputs.  Lines  12  and  13  are 
the  assignments  from  the  stochastic  inputs  to  the  local  C  variables  from  Lines 
8  and  9  of  the  input  model.  Lines  14  and  15  correspond  to  the  local  variables 
assignments  in  Lines  10  and  11  of  the  C  model.  Finally,  Line  16  is  derived  from 
the  ASSERTO  statement  on  Line  13  of  the  C  model. 

Monte-Carlo  Simulation.  The  final  step  of  OSMOSIS  is  the  Monte-Carlo 
simulation  to  estimate  praw  using  (7).  Each  Bernoulli  trial  in  this  simulation 
is  conducted  by  directly  executing  the  dynamically  loadable  executable  of  the 
model.  The  model  source  file  is  compiled  by  gcc,  dynamically  loaded,  then  re¬ 
peatedly  called  for  each  trial.  Before  each  execution  a  random  vector  x  ~  g*  (x) 
is  generated  as  described  in  the  previous  section  and  inititalized.  A  global  flag 
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Name  In  Lmax/G 

p *  i /p* 

dReal  Calls 
none  1  2  1+2 

Time 

none  1+2 

simple  2 

5.859  x  10~3  169 
2.197  x  IIP3  455 

49  38  26  26 
73  57  40  40 

0.15  0.1 

0.21  0.1 

hockey  2 

3.516  x  1CT2  28.4 
1.148  x  ltr2  87.1 

255  213  142  137 

391  328  214  211 

315  228 

364  255 

v  wr  «  10/4 

bcLCKOII  0  12/4 

1.797  x  IIP1  5.6 
1.797  x  KT1  5.6 

479  451  240  240 

1583  1551  792  792 

33  14 

61  28 

V  9  W1 

bounce  2  12/1 

2.997  x  1(T2  33 
1.221  x  lCT2  81 

117  86  59  59 
221  163  111  111 

91  53 

150  84 

Table  1.  AIF  Generation  Results;  In=number  of  inputs;  Lmax=recusrion  depth  limit; 
G=variabfe  grouping  factor,  Time=generation  time  in  seconds;  none,  1,  2  and  1+2 
indicate  which  optimizations  were  used. 


variable  indicating  success/failure  is  also  cleared.  The  function  INPUT_D()  in¬ 
dexes  and  returns  a  value  from  that  vector.  The  ASSERT  ()  statement  tests  the 
condition,  and  if  the  condition  fails  it  sets  the  global  flag  to  true  and  returns. 
Success  or  failure  of  the  trial  is  recorded  based  on  the  value  of  the  flag  variable. 
Trials  resulting  in  an  ASSERT ()  fail  are  inputs  where  IM^$(x)  =  1,  and  those 
where  the  ASSERT ()  does  not  fail  are  inputs  where  Im\=<p(x)  =  0- 

6  Results 

To  evaluate  our  technique,  we  tried  OSMOSIS  on  the  following  example  problems: 
simple  The  example  problem  from  Figure  6a. 

hockey  A  hockey  puck  gets  a  random  initial  impulse  from  a  random  direction. 

Failure  means  that  it  stops  on  a  point  target  after  zero  or  more  bounces, 
backoff  An  exponential  backoff  problem  in  which  two  senders  attempt  up  to  3 
communications.  Failure  occurs  if  transmission  for  either  exceeds  a  deadline, 
bounce  A  ball  is  launched  at  a  random  initial  angle  and  velocity.  We  test  if  it 
falls  in  a  small  hole  after  potentially  bouncing  a  number  of  times. 

Each  of  these  problems  has  the  characteristic  that  the  failure  region  is  disjoint 
in  the  input  space.  For  example  in  the  hockey  problem  there  are  multiple  paths 
by  which  the  puck  can  reach  the  target.  All  experiments  were  peformed  under 
Linux  Ubuntu  12.04  on  a  2.2GHz  Intel  Core  i7  machine  with  16  Gb  of  RAM. 

Table  1  shows  the  results  for  AIF  generation  on  each  example.  For  each 
problem  we  adjust  the  recursion  depth  limit  and  the  variable  grouping  factor 
(number  of  successive  times  each  input  is  split  while  recursing).  We  used  a  larger 
G  for  the  “backoff”  example  because  we  observed  that  a  higher  G  leads  to  better 
performance  for  models  with  larger  numbers  of  inputs.  Recall  from  (12)  that  1/p* 
is  an  estimate  for  -W.  Note  that  the  value  of  p*  value  can  be  used  to  dynamically 
limit  the  recursion  depth  when  generating  the  AIF,  terminating  when  we  have 
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Name  RE  Lmax/G 

p  N  N/N* 

Time  (sec.) 
SMC  total 

CMC 

0.01  10/1 
•ample  12/1 

5.95  x  10~4  1.68  x  107 

5.89  x  10~4  8.95  x  104  187 

6.03  x  10~4  2.64  x  104  636 

6  6 

<0.1  0.1 

<0.1  0.1 

simple  CMC 

0.001  10/1 
12/1 

5.910  x  10-4  1.69  x  10” 

5.910  x  10-4  8.92  x  106  189 

5.910  x  10-4  2.72  x  106  304 

580  580 

4  4.1 

1  1.1 

CMC 

0.01  10/1 
hnrkov  ^ 

6.18  x  10~4  1.58  x  107 

6.18  x  10~4  5.59  x  105  28.3 

6.22  x  10~4  1.74  x  105  90.1 

6.8  6.8 

0.3  228.3 

0.1  255.1 

hockey  CMC 

0.001  10/1 
12/1 

6.215  x  10"4  1.61  x  10” 

6.214  x  10-4  5.56  x  107  29.0 

6.212  x  10-4  1.74  x  107  92.5 

687  687 

25  253 

8  263 

CMC 
0.01  10/4 

bacboff  12/4 

1.21  x  10~4  8.24  x  107 

1.20  x  10~4  1.50  x  107  5.5 

1.21  x  10~4  1.50  x  107  5.5 

25  25 

6  20 

6  34 

backott  CMC 

0.001  10/4 

12/4 

1.193  x  10"4  8.38  x  10” 

1.190  x  10-4  1.51  x  109  5.5 

1.194  x  10-4  1.50  x  109  5.6 

2593  2593 

553  567 

543  571 

CMC 
0.01  10/4 

turner  12/4 

2.96  x  10^b  3.337  x  10s 

3.00  x  10~5  8.464  x  106  39 

2.97  x  10~5  4.104  x  106  81 

133  133 

4.1  57.1 

2.0  86.1 

bounce  CMC 

0.001  10/4 

12/4 

2.989  x  10_b  3.345  x  10iu 

2.993  x  10-5  8.474  x  108  39.5 

2.994  x  10-5  4.068  x  108  82 

13,619  13,619 

432  485 

209  293 

Table  2.  SMC  Results;  RE  =  RE(p)= target  relative  error;  G=group  factor. 


acheived  a  sufficient  gain,  or  when  there  is  insufficient  improvement  from  one 
level  to  the  next.  The  four  columns  under  “dReal  Calls”  show  the  number  of 
calls  that  were  made  to  dReal  using  no  optimization,  using  Optimization  1  only, 
using  Optimization  2  only  and  using  both  optimizations  (see  Section  4.2). 

We  see  that  both  optimizations  are  effective  at  reducing  the  number  of  calls, 
but  that  Optimization  2  peforms  better,  reducing  the  number  of  calls  as  well  as 
total  SMT  solving  time  by  half  in  most  cases.  Also,  while  there  is  some  benift  to 
using  both  optimizations  together,  the  additional  advantage  is  relatively  small. 
This  is  because  when  using  both  optimizations  together,  Optimization  1  can  only 
be  applied  when  the  counter-example  employed  by  Optimization  2  falls  on  a  cube 
boundary,  or  when  analysis  of  a  parent  cube  timed  out  and  is  UNKNOWN. 

Finally,  the  “Time”  column  shows  the  time  to  generate  in  seconds. 

Times  using  no  optimization  (none),  and  using  both  optimizations  (1+2)  are 
shown  to  demonstrate  the  impact  of  the  optimization  techniques.  Note  that  in 
our  current  implementation,  we  do  not  parallelize  the  calls  to  dReal,  which  could 
lead  to  additional  gains. 
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Table  2  shows  the  results  from  the  SMC  phase  of  OSMOSIS.  For  each  sample 
problem,  we  show  the  results  for  target  relative  errors  (RE)  of  0.01  and  0.001. 
At  each  target  RE,  we  compare  Crude  Monte-Carlo  (CMC)  with  SIS  using  two 
different  recursion  depth  limits  as  shown  in  the  Lmax/G  column.  The  probability 
estimate  for  each  experiment  is  shown  in  the  p  column.  We  see  that  the  estimates 
for  CMC  and  SIS  are  very  close  for  each  problem,  and  that  as  expected  the 
agreement  for  those  at  a  relative  error  of  0.001  are  in  closer  agreement. 

The  column  labeled  N  shows  the  number  of  samples  needed  to  acheive  the 
target  relative  error  for  each  experiment,  and  the  column  labeled  N/N*  shows  the 
improvement  of  SIS  over  CMC  for  the  SIS  experiments.  We  can  see  improvements 
ranging  from  a  factor  of  5  to  a  factor  of  over  600.  When  we  compare  the  measured 
N/N*  to  the  values  predicted  by  1/p*  in  Table  1,  we  see  good  agreement.  For 
example,  in  the  “hockey”  problem  with  a  recusion  depth  of  10,  we  got  28.4  as  the 
predicted  improvement,  compared  to  measured  improvements  of  28.3  for  a  target 
RE  of  0.01  and  29.0  for  a  target  RE  of  0.001.  Note  that  the  predictions  for  the 
“simple”  example  slightly  underestimate  the  measured  improvements  because 
the  p*  values  are  close  the  measured  probabilities  while  the  predicition  assumes 
p*  p  and  will  tend  to  underestimate  as  p*  approaches  p. 

That  last  two  columns  show  the  verification  time  for  the  SMC  phase  alone, 
and  for  the  total  time  including  the  abstract  indicator  function  generation  time 
shown  in  Table  1.  We  see  that  SIS  outperforms  CMC  in  all  cases  but  one,  often 
by  multiple  orders  of  magnitude.  Also  since  the  cost  for  generating  the  abstract 
indicator  function  is  fixed  regardless  of  the  target  RE,  there  will  always  be  some 
target  RE  for  which  SIS  outperforms  CMC. 

7  Conclusion 

Statistical  model  checking  (SMC)  is  a  prominent  approach  for  rigorous  analysis 
of  stochastic  systems  using  Monte-Carlo  simulations.  In  this  paper,  we  devel¬ 
oped  a  new  technique,  called  Semantic  Importance  Sampling  (SIS),  to  advance 
the  state-of-the  art  in  applying  SMC  to  compute  the  probability  of  a  rare  event 
using  a  small  number  of  simulations.  SIS  uses  the  semantics  of  the  target  system 
to  recursively  compute  an  abstract  indicator  function  (AIF),  which  is  subse¬ 
quently  employed  to  perforin  SMC.  We  also  present  two  optimizations  to  SIS 
that  reduce  the  number  of  calls  to  SMT  solvers  needed  to  compute  the  AIF. 
We  have  implemented  SIS  in  a  tool  called  OSMOSIS,  and  experimented  with  a 
number  of  examples.  Our  results  indicate  that  SIS  reduces  cost  of  SMC  by  or¬ 
ders  of  magnitude,  and  our  optimizations,  in  combination,  reduce  the  cost  of 
SMT  solving  by  half.  We  believe  that  extending  SIS  to  analyze  stochastic  sys¬ 
tems  compositionally,  and  combining  it  with  symbolic  simulation  techniques,  are 
important  directions  for  future  research. 
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